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ABSTRACT 

We have constructed a grid of non-LTE disk models for a wide range of black hole mass and mass 
accretion rate, for several values of the viscosity parameter a, and for two extreme values of the black 
hole spin: the maximum-rotation Kerr black hole, and the Schwarzschild (non-rotating) black hole. Our 
procedure calculates self-consistently the vertical structure of all disk annuli together with the radiation 
field, without any approximations imposed on the optical thickness of the disk, and without any ad hoc 
approximations to the behavior of the radiation intensity. The total spectrum of a disk is computed 
by summing the spectra of the individual annuli, taking into account the general relativistic transfer 
function. 

The grid covers nine values of the black hole mass between M = 1/8 x 10^ and 32 x 10^ M© with a 
two-fold increase of mass for each subsequent value; and eleven values of the mass accretion rate, each a 
power of 2 times 1 Mq yr~^. The highest value of the accretion rate corresponds to the total luminosity 
L/L-Edd ~ 0.3. We show the vertical structure of individual annuli within the set of accretion disk models, 
along with their local emergent flux, and discuss the internal physical self-consistency of the models. 
We then present the full disk-integrated spectra, and discuss a number of observationally interesting 
properties of the models, such as optical/ultraviolet colors, the behavior of the hydrogen Lyman limit 
region, polarization, and the number of ionizing photons. Our calculations are far from definitive in 
terms of the input physics, but generally we find that our models exhibit rather red optical/UV colors. 
Flux discontinuities in the region of the hydrogen Lyman limit are only present in cool, low luminosity 
models, while hotter models exhibit blueshifted changes in spectral slope. 

Subject headings: accretion, accretion disks — galaxies:active — galaxies muclei 



1. INTRODUCTION 

A black hole that steadily accretes gas from its surround- 
ings at high rates, but below the Eddington limit, should 
form a geometrically thin accretion disk supported by the 
residual angular momentum of the gas. For black hole 
masses in the range 10^ — 10^ Mq, thought to be typ- 
ical of bright active galactic nuclei (AGN) and quasars, 
the peak effective temperature is expected to be around 
10'' — 10^ K. This estimate is roughly consistent with the 
observation that quasar spectra peak in the ultraviolet, 
bolstering the belief that quasars are indeed powered by 
accretion onto massive black holes. 

Over the years, many authors have attempted to cal- 
culate detailed theoretical spectra of geometrically thin, 
optically thick accretion disks in order to compare with 
observation, with varying degrees of sophistication (e.g. 
Kolykhalov & Sunyaev 1984; Sun & Malkan 1989; Laor & 
Netzer 1989; Ross, Fabian, & Mineshige 1992; Shimura & 
Takahara 1993, 1995; Dorrer et al. 1996; and SinceU & 
Krolik 1998). As reviewed recently by Krolik (1999a) and 
Koratkar & Blaes (1999), these and other models gener- 



ally suffer a number of problems when trying to simulta- 
neously explain various features of the observations, e.g. 
optical/ultraviolet continuum spectral shapes, the lack of 
observed features at the Lyman limit of hydrogen, polar- 
ization, the origin of extreme ultraviolet and X-ray emis- 
sion, and correlated broadband variability. It may be that 
the resolution of these problems requires drastic revision of 
the accretion disk paradigm. Alternatively, it may simply 
be that the theoretical modeling to date is still too crude 
to do justice to the inherent complexities of the accre- 
tion flow. In particular, inclusion of non-LTE effects and 
detailed opacity sources, Comptonization, and interaction 
between the disk and X-ray producing regions should all 
be taken into account. 

We have embarked on a long term program to construct 
detailed model spectra of accretion disks in the axisym- 
metric, time-steady, thin-disk approximation. At the peak 
temperatures, the most important opacity is provided by 
electron scattering, but bound-free and free-free contin- 
uum opacities due to hydrogen and helium can also be sig- 
nificant, so we include these opacities in this study. Since 
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scattering can be the dominant opacity and the densities 
can be rather low, departures from local thermodynamic 
equilibrium (LTE) can be significant, so we include non- 
LTE effects as well. We include the effects of relativity on 
the disk structure and on the transport of radiation from 
the disk to infinity. We construct a grid of models for a 
black hole spin of a/M = (Schwarzschild black hole), and 
0.998 (maximum rotation Kerr black hole), luminosities 
between 3 x 10~^ -^Edd and 0.3 1/Edd, and black hole masses 
between Mg = IQ-^M/Mq = 0.125 to Afg = 32. To de- 
termine the surface mass density of the disk, we assume 
that the viscous stress scales as tr4, = aPtotai (Shakura & 
Sunyaev 1973), where tr<t, and Ptotai are the vertically inte- 
grated viscous stress and total vertically-integrated pres- 
sure, respectively. 

In previous papers (Hubeny & Hubeny 1997 and 1998a, 
hereafter Papers I and II respectively), we presented de- 
tailed spectra and the vertical structure of individual an- 
nuli. In this paper we integrate such spectra over radius 
for a grid of 99 mass/accretion rate combinations appro- 
priate for quasars. We do not include annuli with low 
effective temperatures (Tes < 4000 K) as these require 
molecular opacities for accurate computation. The spec- 
tra of hot annuli can be affected by Compton scattering, 
which we have not included in the calculation, but will 
in future work. It is likely that metal opacities will also 
modify the final spectrum (cf. Hubeny & Hubeny 1998b), 
so wc consider this work as a benchmark for future metal 
line-blanketed models. 

This paper is organized as follows. In section 2 we de- 
scribe our model assumptions and computational methods. 
Then in section 3, which constitutes the bulk of the paper, 
we present our results. We start by showing the vertical 
structure of individual annuli within the set of accretion 
disk models, along with their local emergent flux. We then 
discuss the internal physical self-consistency of these mod- 
els, before presenting the full disk-integrated spectra. We 
finish section 3 with a discussion of a number of observa- 
tionally driven issues: optical/ultraviolet colors, spectra in 
the hydrogen Lyman limit region, polarization, and ioniz- 
ing continua. Finally, in section 4 we summarize our con- 
clusions, in particular pointing out the additional physics 
which will be included in future papers of this series. 

2. MODEL ASSUMPTIONS AND COMPUTATIONS 

To construct a model of an accretion disk, we assume 
the vertical disk structure can be well approximated by 
one-dimensional equations; that is, we assume the disk is 
locally plane parallel. We assume that, on average, the 
disk is static in the corotating frame (in reality disks are 
subject to many instabilities which invalidate this approx- 
imation), and that the only energy transport is due to 
radiation flux in the vertical direction, i.e. we ignore con- 
vection and conduction. By assuming time-steadiness and 
local radiation of dissipated heat, we can write down and 
solve the equations for the disk structure (Page & Thorne 
1974): these equations are summarized in Paper II. For a 
given radius r, the (one-sided) flux (and thus the effective 
temperature, Tes) is determined by 



F{r) = ubT; 
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where (Tb is the Stefan-Boltzmann constant, M the mass 
accretion rate, and is a relativistic correction factor 
(Page & Thorne 1974, in the notation of Krolik 1999a). 
Following the usual practice for geometrically thin accre- 
tion disks, we assume that there is no torque at the inner- 
most stable circular orbit in all our disk models. There are 
reasons to question this assumption (Krolik 1999b, Gam- 
mie 1999); if it fails, the disk spectrum and polarization 
could be substantially changed (Agol & Krolik 1999). 

We have calculated models for two values of the viscos- 
ity parameter: a — 0.01, and a = 0.1. The choice of 
a = 0.01 is near the value expected from simulations of 
the magneto-rotational instability in accretion disks (e.g. 
Balbus fc Hawley 1998), while a = 0.1 represents a typi- 
cal value of the viscosity parameter used in other studies. 
Smaller a or a stress which scales in proportion to the 
gas pressure would lead to a larger surface mass density, 
higher density, and, presumably, spectra closer to LTE. 

If the disk rotates on cylinders, and the angular fre- 
quency of each cylinder is the one appropriate to a circular 
orbit at the midplane, there is a vertical component to the 
effective gravity proportional to height z above the disk 
midplane. We ignore the self-gravity of the disk, an excel- 
lent approximation for the radii important to our problem. 

We have assumed that the local dissipation is propor- 
tional to the local density, except for the top 1% of the disk 
where we force the dissipation to decline - see Paper II. 
This distribution is chosen in order to yield a hydrostatic 
equilibrium solution in the bulk of the disk when radiation 
dominates (cf. Shakura & Sunyaev 1973) while retaining 
the possibility of thermal balance in the outer-most lay- 
ers even in the absence of Comptonization. In Paper II 
(in particular, see figures 11 and 12 there) we showed that 
neither the choice of the division point between the regions 
of vertically constant and declining viscosity, nor the slope 
of the power law for the viscosity in the declining regime, 
changes the predicted continuum emergent spectrum sig- 
nificantly. Disks with this structure can be convectively 
unstable, however. Convection can be expected to accel- 
erate heat loss, leading to a disk structure that is rather 
thinner and denser (Bisnovatyi-Kogan & Blinnikov 1977), 
although perhaps not substantially so (Shakura, Sunyaev 
& Zilitinkevich 1978). Even annuli with entropy increas- 
ing upward, so that they are stable according to the usual 
Schwarzschild criterion, may still nevertheless be subject 
to convectivc instabilities mediated by thermal conduction 
along weak magnetic field lines (Balbus 1999). How this 
would manifest itself in the presence of turbulence driving 
the radial angular momentum transport is not clear, how- 
ever. In any case, as stated above, we completely ignore 
convectivc heat transport in all our models here. 

When either radiation or gas pressure dominates, the 
surface mass density may be readily calculated in the grey, 
one-zone, diffusion approximation. However, when the two 
are comparable to each other, the disk structure equations 
combine to form a single tenth-order polynomial equation 
in one of the variables (e.g., S^/^). We solve this equation 
iteratively, using the two limiting regimes to bracket the 
root. 

Once the effective temperature, surface mass density, 
and dissipation profile are determined, the vertical disk 
structure can be computed. To compute the vertical struc- 
ture of a given annulus, we solve simultaneously the entire 
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set of structural equations: hydrostatic equilibrium in the 
vertical direction; local energy balance; radiative transfer; 
and, since we are not generally assuming LTE, statisti- 
cal equilibrium for all selected energy levels of all selected 
atoms and ions. The equations are solved for the whole ex- 
tent of the disk between the midplane and the surface using 
proper boundary conditions. For details, the reader is re- 
ferred to Paper 11. We stress that no ad hoc assumptions 
about the nature of the radiation field or the radiative 
transfer are made; for instance, we do not use the diffu- 
sion approximation or an escape probability treatment, or 
any assumption about the angular dependence of the spe- 
cific intensity; the radiative transfer is solved exactly. To 
represent the radiation field, we use about 150 frequency 
points placed to define all continuum edges and resolve 
any frequency-dependent structure in the emergent inten- 
sity. The minimum frequency is set to 10^^ Hz, while the 
highest frequency is chosen so that even at the midplane 
the intensity for frequencies higher than the maximum is 
negligible. In terms of the midplane temperature (Paper 
II and Hubeny 1990) 

T^nid « Teff(3r^id/8)^/^ = Teff[(3/8)(S/2)(xii)]^/^ , (2) 

this goal can be achieved by setting the maximum fre- 
quency to 

i^max = 17A:Tmid/^- (3) 

Here XR is the Rosseland mean opacity (per gram) , which 
we take for simplicity to be 0.34, the value corresponding 
to an opacity dominated by electron scattering in a fully 
ionized H-He plasma of solar abundance. 

We consider disks composed of hydrogen and helium 
only. We include metals in computing the molecular 
weight, but for this work, we ignore their effects on the 
opacity. Hydrogen is represented essentially exactly: the 
first 8 principal quantum numbers are treated separately, 
while the upper levels are merged into a single non- 
LTE level accounting for level dissolution as described by 
Hubeny, Hummer, & Lanz (1994). We do, however, as- 
sume complete Z-mixing; given the high electron density 
in these environments, this should be a good approxima- 
tion. Neutral helium is represented by a 14-level model 
atom, which incorporates all singlet and triplet levels up 
to n = 8. The 5 lowest levels are included individually; sin- 
glet and triplet levels are grouped separately from n = 3 to 
n = 5, and we have formed three superlevels for n = 6, 7, 
and 8. The first 14 levels of He+ are explicitly treated. We 
assume a solar helium abundance, A''(He)/A'^(H) = 0.1. 

The opacity sources we include are all bound-free tran- 
sitions (continua) from all explicit levels of H, He I, and 
He II; free-free transitions for all three ions, and electron 
scattering. For the coolest models (Teg < 9000 K), we 
also consider the H~ bound-free and free-free opacity, as- 
suming LTE for the H^ number density. In this paper, 
we assume coherent (Thompson) scattering. As was dis- 
cussed in Paper II, effects of non-coherent (Compton) scat- 
tering are negligible for models with T^ff around or below 
10^ K. Therefore, most models of the present grid (see 
Sect. 3) are not influenced by the effects of Comptoniza- 
tion, although the hottest ones may be. We have recently 
implemented Comptonization in our modeling code, and 
checked that this is indeed the case; however, we choose 



to neglect Comptonization in this paper in order to pro- 
vide a benchmark grid of models computed using classical 
approximations of H-He composition and without Comp- 
tonization. The effects of Compton scattering will be in- 
cluded in a future paper, where we will also discuss in 
detail its influence on the emergent spectra. It can be 
expected to become especially important when the dissi- 
pation per unit mass is enhanced near the disk's surface. 

The structure equations are highly nonlinear, but are 
very similar to corresponding equations for model stellar 
atmospheres. We use here the computer program TLUS- 
DISK, which is a derivative of the stellar atmosphere pro- 
gram TLUSTY (Hubeny 1988). The program is based on 
the hybrid complete-linearization/accelerated lambda it- 
eration (CL/ALI) method (Hubeny & Lanz 1995). The 
method resembles traditional complete linearization, how- 
ever the radiation intensity in most (but not necessarily 
all) frequencies is not linearized; instead it is treated via 
the ALI scheme (for a review of the ALI method, see 
e.g., Hubeny 1992). Moreover, we use Ng acceleration 
and the Kantorovich scheme (Hubeny & Lanz 1992) to 
speed the solution. We start with a grey atmosphere so- 
lution (Hubeny 1990), flrst solving for the disk structure 
assuming that the statistical equilibrium is described by 
LTE. Using this as a starting point, we drop the LTE as- 
sumption and compute the disk structure using the full 
statistical equilibrium equations, as described in Paper II, 
but assuming that the line transitions are in detailed ra- 
diative balance. In other words, the statistical equilibrium 
equations explicitly contain the coUisional rates in all tran- 
sitions, and radiative rates only in the continuum transi- 
tions. We have considered 70 discretized depth points. 
The top point is set to mi = 10""^ g cm~^, where m is 
the column mass, i.e., the mass in a column above a given 
height. The last depth point is the column mass corre- 
sponding to the midplane, and is given by S/2. The depth 
points are equally spaced in logarithm between these two 
values. The models were computed on a DEC Alpha with 
500 MHz clock speed; with the above values for the num- 
ber of frequency and depth points, the LTE models for 
individual rings required typically 5-10 iterations with ap- 
proximately 2.5 seconds per iteration, while the non-LTE 
models required typically 5-10 (for hotter annuli), or 10- 
30 iterations (for cooler annuli), with about 6 seconds per 
iteration. 

We do not compute here models treating radiative rates 
in line transitions explicitly. We have considered such 
models in Paper II and found that including lines explic- 
itly does not change the vertical structure or emergent 
continuum radiation signiflcantly. These models are also 
much more time consuming. However, the most impor- 
tant physical point is that we have found in test calcula- 
tions that line profiles are influenced significantly by the 
effects of Compton scattering. Since we are neglecting the 
Comptonization here to provide a benchmark grid of clas- 
sical H-He models without Comptonization, we feel that 
including lines at this stage would not be much more than 
a numerical exercise. We therefore defer treatment of more 
realistic models including lines and Comptonization to a 
future paper. 

In the course of calculating atmosphere models, we 
sometimes ran into difficulties with convergence. For very 
low and very high temperatures, we could not get conver- 
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gence, a problem which may be amehorated in the future 
by including more sources of opacity (such as metals at 
high temperatures and molecules at low temperatures). At 
certain radii and certain depths within the disk, we found 
ionization fronts in helium around Toff ^ 35, OOOK which 
are very narrow in extent and very sensitive to the temper- 
ature within the disk, causing limit cycle behavior where 
the helium alternates between being mostly doubly ionized 
or recombined to a singly ionized stage during successive 
steps, preventing solution of the atmosphere structure. We 
solved this problem by computing the disk structure for 
radii just smaller or larger than the radii where the front 
exists, and then using these solutions as starting solutions 
for the radii at which the fronts exist. In practice, the 
range of radii where this problem exists is narrow so that 
the uncertainty in the structure does not affect the overall 
disk spectrum. We also found similar He I/He II ioniza- 
tion fronts at lower temperatures, Tcff ~ 15, OOOK, and for 
hydrogen around and below 9000 K. 

To find the total disk spectrum, we divided the disk into 
25-35 radial rings, spaced (roughly) logarithmically. At 
each ring, after computing the vertical structure, we per- 
form a detailed radiation transfer solution for the Stokes 
vector as a function of frequency and angle. The spec- 
trum is found by integrating the total emergent intensity 
over the disk surface using our relativistic transfer function 
code (Agol 1997). The transfer function computes the tra- 
jectories of photons from infinity to the disk plane, finding 
the emitted radius, redshift, and intensity at each image 
position at infinity for a given observation (Cunningham 
1975). In this paper, we neglect the effects of radiation 
which returns to the accretion disk. 



3. RESULTS 

The parameter space of our grid of models is displayed in 
figure H The defining parameters are Mg, the black hole 
mass expressed in IO^AIq, and M, the accretion rate in 
units of Mq yi:~^. The grid covers nine values of the black 
hole mass between Mg = 1/8 and 32; each subsequent 
mass is twice the previous mass. An analogous approach is 
used for the mass accretion rate, i.e. eleven values for each 
black hole mass which are powers of 2 times 1 Mq yr~^. 
The highest value of the accretion rate in each case is cho- 
sen to make L/L^dd = 0.286; i.e., M{Mqji-'^) = 2Mq. 
In the following text, we refer to this highest value as 
L/ L^dd = 0.3. the ten subsequent values have half the ac- 
cretion rate (and luminosity) of the previous model. Our 
grid spans a range similar to models which have been pre- 
viously used to fit quasar spectra (Sun & Malkan 1989; 
Laor 1990). 

Our basic grid assumes a/M — 0.998. However, we 
have also constructed a parallel grid for a Schwarzschild 
black hole with the same black hole masses and values 
of L/L-Edd- The mass accretion rates for a/M = are 
a factor of 5.613 higher than the corresponding values for 
a/M = 0.998 because the radiative efficiency of a Kerr disk 
is that much greater (e.g. Shapiro & Teukolsky 1983). 
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Fig. 1. — The parameters of our grid of disk models for the Kerr 
hole case. Upper panel: -L/-LEdd versus black hole mass, expressed 
in 10® Mq. Middle panel: the mass accretion rate, M, expressed in 
AIq yr^-"^. Lower panel: the maximum effective temperature (which 
is found at r ~ 1.5 rg). 

For all disk models, we compute detailed vertical struc- 
ture at the following radii (expressed as r/rg, where rg = 
GM/c^ is the gravitational radius): r/rg = 1.5, 1.7, 2, 2.5, 
3, 3.5, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 50, 
60, 70, 80, 90, 100, 120, 140, 160, 180, 200, 250, 300, 400, 
500; that is, while the corresponding effective tempera- 
ture Tcff > 4000 K. For the Schwarzschild grid, we use the 
above values of r/r^, multiplied by 5. The actual number 
of computed annuli thus depends on the basic parameters 
of the disk. 

3.1. Properties of Models 

In this section, we discuss the behavior of individual an- 
nuli, while in the rest of the paper we will present the emer- 
gent radiation integrated over the whole disk. We have 
chosen a model with Mg = 1, M = 1 (i.e., L/L^dd = 0.15) 
as a representative model for displaying various quanti- 
ties. The behavior of other individual disk models is sim- 
ilar. Figure || displays the local electron temperature as a 
function of position for all individual annuli. The position 
is expressed as column mass, m, above the given depth. 
As mentioned in Sect. 2, the uppermost point was chosen 
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to be m = 10^'^ g cm^^, while the highest value corre- 
sponds to the midplane of the disk. The temperature is 
a nearly monotonic function of depth, although there is a 
slight temperature rise at the surface for some models. The 
detailed behavior of temperature for several representative 
annuli was discussed extensively in Paper IL The behavior 
of temperature for all annuli is easily understood. In the 
LTE approximation, the midplane temperature is given 
by equation (Q), i.e. it is proportional to Tcs{r)[mo{r)]^^'^ , 
where mg = S/2 is the total column density at the mid- 
plane. Neglecting for simplicity the relativistic corrections. 
Toff oc r"'^/'' - see equation (||), while the column density 
mg cx r'^/^ for radiation pressure dominated annuli (see 
Eq. 18 of Paper II), which is the case for the models 
considered in figure y. Therefore, the LTE approximation 
predicts that the midplane temperature, Tmid oc r^"^/®. In 
contrast, the surface temperature is proportional to Tofr; 
therefore, Tgurf oc r~^/**. Figure || shows that these scal- 
ings do in fact hold approximately. The ratio of the lowest 
and highest radii is roughly 45, so if the LTE approxima- 
tion held, the surface temperature would vary by a factor 
of 17 at the surface, and by a factor of 4.2 at the midplane. 
The model values are 14 and 4.6, respectively. 
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Fig. 2. — Temperature as a function of depth for the individual 
annuli. The curves correspond to radii r/rg (from top to bottom) 
1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 7, 8, 9, 10, 12, 14, 16, 18, 20, 25, 30, 40, 
50, 60, 70, 80, and 90. 

Figure ^ displays the run of mass density for the individ- 
ual annuli. The central density is lowest for the inner an- 
nuli, and increases with increasing distance from the black 
hole, while the density close to the surface exhibits a more 
or less reverse behavior. Note that the disk becomes opti- 
cally thin below a column mass of m ~ Xr^ — 2-9 g cm~^, 
where the density is substantially lower (in some cases by 
orders of magnitude) than the midplane density. This 
should be borne in mind when comparing our results to 
those of previous workers who assumed constant density 
slabs (e.g. Laor & Netzer 1989). An explanation of the 
behavior of the density is given in the Appendix. 
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Fig. 3. — Mass density as a function of jdepth for the individual 
annuli. The radii are the same as in Fig. H. At the midplane, the 
highest density corresponds to the largest r/vg (equal to 90), while 
the lowest density corresponds to r/rg = 2.5. The two uppermost 
curves for m < 10"^ correspond to r/rg = 1.5 and 2, respectively. 

A more interesting behavior is exhibited by the H I 
and He II ground state number densities, displayed in fig- 
ures ^ and 1^. For hot, inner annuli, hydrogen remains 
ionized throughout the entire vertical extent of the disk, 
while starting with the annulus at r = 40rg (with T^s « 
15, OOOK), there is an appreciable portion of neutral hydro- 
gen in the outer layers (m < 10^ g cm~^). Consequently, 
the Lyman edge opacity for these latter models becomes 
rather large (see below) . Similarly, helium exhibits a tran- 
sition, between r = I2rg (with T^s ~ 36, OOOK) and 
r = lirg (with T^g ~ 32, 500K) from being predominantly 
doubly ionized, to a situation where the single ionized he- 
lium is a dominant stage of ionization at some outer part of 
the disk. For the annuh with r < 40rg (with Tcs > 15, 000 
K), helium remains singly ionized throughout the entire 
upper part of the disk, while for more distant, cooler an- 
nuli, neutral helium becomes more and more dominant in 
the outer layers. 



M = 1; Mdot = 1; L/L = 0.15 
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Fig. 4. — Number density of the ground state of H I as a function 
of depth (full lines), and the proton number density (dotted lines). 
The radii are the same as in Fig. bl Selected models are labeled by 
the radial coordinate, r/rg, and the models for r/vg = 30 and 40, 
which delimit the regions between the full and partial ionization of 
hydrogen, are drawn by bold lines. 
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Fig. 5. — Number density of the ground state of He II as a function 
of depth. The radii are the same as in Fig. U. Selected models are 
labeled by the radial coordinate, r/vg, and the models for r/vg = 12, 
14, and 40, are drawn by bold lines. 

Figure ^ displays the emergent flux for all annuli. The 
upper panel shows the non-LTE models, while the lower 
panel shows the LTE models. The behavior of the emer- 
gent flux is analogous to that discussed at length in Papers 
I and II. Non-LTE models exhibit the He II Lyman edge 
in emission for r <\2rg (with T^s > 36,000K), while for 
more distant annuli the flux in the He II Lyman contin- 
uum drops to almost zero. This is a consequence of the 
transition from the dominance of double-ionized helium 
to singly-ionized helium, as displayed in figure |^. Anal- 
ogously, the neutral- hydrogen annuli for r > 40 r^, with 
Teff > 15, 000 K, (see figure 0) show a significant hydrogen 
Lyman edge in absorption. 

A comparison between LTE and non-LTE results reveals 
several interesting effects: in LTE, the He II Lyman edge 
appears to be in weaker emission or even in absorption 
for hot, doubly-ionized, annuli. Most importantly, non- 
LTE effects reduce the hydrogen Lyman edge (cf. Sun 
& Malkan 1989, Shields & Coleman 1994, Storzer et al. 
1994, Papers I and II): the edge predicted in non-LTE 
models is typically in weaker emission for emission edges, 
and in weaker absorption for absorption edges. We dis- 
cuss the behavior of the hydrogen Lyman edge in the full 
disk-integrated spectra in more detail in Sect. 3.5 below. 

Finally, figure M shows a comparison of predicted local 
non-LTE flux ana the blackbody flux corresponding to the 
same effective temperature. As already discussed in Pa- 
per I, the non-LTE models exhibit a much wider spectral 
energy distribution than blackbodies. Consequently, an 
often-used mapping of a given spectral region to a certain 
radial position in the disk, which may be used for a black- 
body flux because of its sharp variation with frequency, is 
much less satisfactory for more realistic non-LTE vertical 
structure models. Another crucial feature of the present 
models is that compared to blackbodies the energy is re- 
distributed. At low temperatures, absorptive opacity is 



relatively important, so that flux is shifted from frequen- 
M = 1; Mdot = 1; L/L = 0.15; non-LTE 




Fig. 6. — Predicted local emergent flux (in erg cm Hz~^) for 
the individual annuli for non-LTE models (upper panel), and LTE 
models (lower panel). The models for r/rg = 12, 30, and 40, are 
drawn by bold lines. 




frequency 



Fig. 7. — Predicted local emergent flux (in erg cm Hz~^) for 
the individual annuli for non-LTE models (full lines), together with 
a blackbody flux corresponding to the same effective temperature 
as the given annulus (dotted lines). 

cies where the opacity is high to frequencies where it is low; 
i.e., away from ionization edges. Another way of viewing 
this effect is to note that in bands where the opacity is 



7 



low, escaping photons are created deeper inside the disk 
where the temperature is higher. At high temperatures, 
the opacity becomes scattering dominated. In these rings, 
there is a trade-off between the effect just mentioned (lower 
opacity bands allow one to see photons created at higher 
temperature) and scattering blanketing (which retards the 
escape of photons created deep inside). The latter effect 
is the one that creates a "modified blackbody" spectrum 
when there is no temperature gradient. For much of the 
parameter range of interest, the net effect is to shift flux 
toward higher frequencies. 

3.2. Consistency of Models 

Having computed the detailed vertical structure of the 
annuli, we have to address the question of self-consistency 
of the model assumptions. 

First, we check that the computed disks are indeed ge- 
ometrically thin, i.e., the disk height, H_, is much smaller 
than the radial coordinate r. In figure M we show the ra- 
tio of the disk height to the radial coordinate, H/r, as a 
function of the radial coordinate (in units of gravitational 
radius), for a disk with Mg = 1, and for various values of 
M (or luminosity) . The behavior of disks for other values 
of the black hole mass is almost identical. (This is expected 
for radiation pressure supported, electron scattering dom- 
inated disks, because then H/r can be written as L/ L^^^ 
times a function oir /rg, cf. equation 53 of Paper II.) Only 
for the most luminous disk does the ratio H/r approach 
10% at r/rg — 2.5. The height of other annuli is lower, 
and the maximum height for less luminous disks is progres- 
sively lower. These results are in good agreement with the 
older models of Laor & Netzer (1989, cf. their figure 1 and 
the surrounding discussion). We therefore conclude that 
our disks do indeed satisfy the thin disk approximation. 



M =1; L/L 



0.3 (top) - 0.002 (decreasing by 2) 
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Fig. 8. — The ratio of the disk height to the radial coordinate, 
H/r, as a function of the radial coordinate (in units of gravitational 
radius), for a disk with Mg = 1, and for various values of M (or 
luminosity). The crosses represent the actual computed models of 
the individual annuli. 

Another important concern is the presence of vertical 
density inversions within the disk. Since we neglect con- 
vection, sharp temperature gradients and density inver- 
sions are found at ionization transitions occurring in re- 
gions where gas pressure contributes significantly to sup- 



port against gravity. In figure || we display the structure of 
several annuli that illustrate this behavior. The hottest an- 
nulus shown in that figure {r/rg = 50) has Toff — 9400 K; 
it is the outermost ring with no density inversion. Cooler, 
more distant annuli show a progressively stronger inver- 
sion. The inversion is created by the abrupt fall in electron 
density when H recombines; in order to maintain the pres- 
sure gradient required for hydrostatic balance, the mass 
density must increase to compensate. A sharp rise of local 
temperature, displayed in the upper panel of figure ^ is 
caused by a rapid increase of opacity when going inward 
due to an ionization front. The temperature is essentially 
a function of the Rosseland mean opacity, T cx T^ogg, but 
since tross is a sharp function of the column mass m, the 
function T{m) also exhibits a sharp increase with m. 



Mdot = 1; L/L 




10" lO"- 10^ 

column mass (g cm~^) 



Fig. 9. — Upper panel: temperature as a function of position for 
a disk with Mg = 2 and M = IMq yr~^, for "cool" annuli (for radii 
50, 60, 70, 80, 90, 100, 120, and 160 rg. The corresponding Teff's 
are 9400, 8300, 7400, 6700, 6200, 5700, 5000, and 4100 K. The lower 
panel: mass density for the same annuli. 

These models should be taken with caution. In the ab- 
sence of detailed hydrodynamical calculations which would 
allow for convective instability and determine the vertical 
structure properly, we do not know what is the correct 
emergent radiation. We can nevertheless obtain a rough 
indication of model uncertainties by comparing the pre- 
dicted flux from the non-LTE models that neglect con- 
vection, and a blackbody flux corresponding to the same 
effective temperature. 

We present in figureJlOl emergent flux for three selected 



annuli shown in figure BTfor r/rg = 50, 90, 160, with corre- 
sponding effective temperature equal to 9400, 6200, 4100 
K, respectively. The predicted flux for non-LTE models 
is shown together with the corresponding blackbody flux. 
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We see that the models at the hotter end of the density in- 
version sequence exhibit an appreciable Balmer edge, while 
the cooler models are reasonably well approximated by 
blackbodies. We thus feel that approximating even cooler 
models, which we do not compute here, by blackbodies is 
probably not worse than other approximations that un- 
derlie the entire calculation. 




frequency 



Fig. 10. — Predicted local emergent flupt. (in erg cm s^^ 
Hz~^) for three selected annuli shown in Fig. bI Full lines: non- 
LTE models; dotted lines: blackbody flux corresponding to the 
same effective temperature as the given annulus. The curves are 
labeled by the value o{ r/vg = 50, 90, 160; the corresponding 
Tcft = 9400, 6200, 4100 K, respectively. 



3.3. Disk-integrated Spectra 

In this section, we present disk-integrated spectra for 
selected disk models. The full grid of models is not pre- 
sented here, but is available to interested researchers upon 
request. 

Note that in all the spectral energy distributions shown 
in this paper, the quantity is the specific luminosity 
that an observer along a particular viewing angle would 
infer the source to have if it were isotropic, i.e. ii Fy is the 
measured specific fiux and d is the distance to the source, 
then = Aird^F^. In the subsequent plots, we display, as 
is customary, the quantity i>L,y, which represents a lumi- 
nosity per unit logarithmic interval of frequency (photon 
energy) . 

As discussed above, we consider the spectra of the outer 
annuli which were not computed (for Toff < 4000 K) to be 
given by the black-body energy distribution corresponding 
to the effective temperature of the annulus. An impor- 
tant parameter is the outer cutoff of the disk. In order to 
avoid problems with an improper choice of the outer cutoff, 
we have chosen the cutoff radius Tout in such a way that 
Teff{rout) is equal to a specific limiting temperature, Tn^. 
We have chosen Tum = 1000 K, which guaranties that the 
total emergent flux at v = 10^*, which is the lowest fre- 
quency considered in our integrated disk spectra, is not 
significantly influenced by cooler annuli with Toff < Tiim 
to within a few per cent. 



M =2; Mdot = 1; L/L =0.07; cos i = 0.8 




frequency (Hz) 



Fig. 11. — Integrated spectral energy distribution uL^, for a disk 
with Mg = 2 and M = I Mq yr"! (L/I/Edd = 0.07), with inclina- 
tion cosi = 0.8, and with various choices of Tout. Solid lines: the 
automatic choice of rout, corresponding to = 1000 K (in this 

case, Tout/rg ~ 1060). Dashed lines: rout/^g = 160; Dotted line: 
rout/rg = 60; the curves are labeled by the value of rout/^g). Bold 
lines correspond to spectra computed for non-LTE models down to 
Tcff = 4000 K; thin lines show the spectra computed by replacing 
the local spectra of annuli with Tf,g < 9000 K by blackbodies. 

Figure O shows the effect of the outer cutoff, as well 
as of the degree of approximation in the treatment of cool 
annuli, for the disk model displayed in figure || (Mg = 2; 
M = IMq yr-\ i.e., L/L^dd = 0.07). Three cutoff radii 
are shown: rout/'^g = 60, which is the radius where the 
density inversion sets in; rout/^'g — 160 (the radius where 
Toff reaches the our minimum value of 4000 K), and fi- 
nally the default value which corresponds to T\i^ — 1000 
K, which in this case happens at rout/rg « 1060. Neglect- 
ing all annuli cooler than 9000 K (dotted line) does not 
influence the flux blueward of the Balmer limit, but the 
optical and IR flux are seriously underpredicted. Neglect- 
ing all non-computed annuli cooler than 4000 K (dashed 
lines) produces the correct flux in the region of the Balmer 
edge, and very nearly the correct flux in the optical range 
{v > Ax lO^** Hz). The effect of uncertainties in the mod- 
els with density inversion (Tefi < 9000 K) is estimated by 
comparing bold lines, which show spectra computed for 
non-LTE models down to Toff = 4000 K, with the corre- 
sponding thin lines, showing the predicted integrated spec- 
tra when all annuli with Toff < 9000 K are assumed to emit 
locally as blackbodies. The maximum difference in the 
integrated spectra is found redward of the Balmer edge; 
the non-LTE models of cool annuli produce an integrated 
flux that is about 16 % higher than that corresponding to 
replacing the Tcff < 9000 K annuli by blackbodies. There- 
fore, the effect is not very large, which gives us confidence 
that our overall integrated spectra are not significantly in- 
fluenced by the uncertainties associated with the density- 
inversions present in the cool annuli. However, the pre- 
dicted feature at the Balmer edge should be viewed with 
caution. 

We first present integrated spectra for a disk of given 
mass and luminosity, Mg = 1, L/L^dd = 0.15 (i.e. 
M — IMq yr~^)i and for various values of the incli- 
nation angle i, ranging from cosi — 0.99 (i.e., the disk 
seen almost face-on), to cosi = 0.01 (the disk seen almost 
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edge-on) - figure We can clearly see the impact of rel- 
ativistic boosting and aberration: the flux at the highest 
frequencies (radiated where the disk is most relativistic) 
is boosted strongly for viewing angles near the disk plane. 
Also for such viewing angles, gravitational light bending 
counteracts the Newtonian "cos0" projection effect. 

M =1; Mdot = 1; cos i = 0.99, 0.8, 0.6, 0.4, 0.2, 0.01 



2 ' - 2^ (increasing by 2); Mdot = 1; cos i = 0.1 




frequency (Hz) 

Fig. 12. — Integrated spectral energy distribution uLi, (in erg 
s~^) for a disk with Mg = 1, M = IMq yr~^, and for various 
inclinations. 

M.= l; Mdot = Z'" - 2' (increasing by 2); cos i = 0.8 




10 10 
frequency (Hz) 



Fig. 13. — Integrated spectral energy distribution uL,j (in erg 
s~^) for a disk with Mg = 1, cosi = 0.8, and for various values of 
the mass accretion rate (luminosity), ranging from L/Z/Edd = 0.3 
(i.e., M = 2Mq yr-l) to L/^Edd = 3 X lO"" (i.e., M = 
yr~^). Full lines: non-LTE models; dotted lines: LTE models. 

In the following, we present the integrated spectra for 
one value of inclination. We chose cosz = 0.8 (i.e., 
i « 37°), which is a value relatively close to face-on, and 
which thus may serve as a typical value for type 1 AGN and 
QSO's based on unification arguments (e.g., Krolik 1999a). 
Figure 13 shows the integrated spectra for one particular 
value oithe black hole mass, Mg = 1, and for various lumi- 
nosities (mass accretion rates) . Full lines display non-LTE 
models, while the dotted lines display LTE model predic- 
tions. The spectral energy distribution is hardest for the 
largest luminosity. The non-LTE effects are most impor- 
tant in the He II Lyman continuum (i/ > 1.36 x 10^^ Hz), 
and also for predicting the detailed shape of the hydrogen 
Lyman edge for intermediate and low luminosities. 



non-LTE (full); LTE (dot) 

33 




frequency (Hz) 

Fig. 14. — Integrated spectral energy distribution uLi, (in erg 
s~^) for a disk with M = IMq yr"-*^, cosi = 0.8, and for various 
values of the black hole mass, Mq, ranging from Afg = 32 to 0.5. 

In figure |lj we display the sequence of predicted spec- 
tra for models with a fixed mass accretion rate (i.e., total 
luminosity) and varying black hole mass. The energy dis- 
tribution is progressively shifted towards more energetic 
photons for lower masses, because disks around less mas- 
sive holes have smaller radiating surface areas. The non- 
LTE effects are important for all disks; for hotter ones the 
largest departures from LTE arc seen in the He II Lyman 
continuum, and in cooler ones in the hydrogen Lyman con- 
tinuum. 

(M , Mdot) pairs with the same T (r/r ); cos i = 0.8 
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Fig. 15. — Integrated spectral energy distribution uLi, (in erg 
s~^) for the pairs of (Mg, M) which have the same effective tem- 
perature distribution, Tcff (r/r^), and which should therefore have a 
rather similar spectrum shape. The curves are labeled by the values 
of (Mg, M). Full Unes: non-LTE models; dotted lines: LTE models. 

Figure ^s] shows the spectral energy distribution for a 
sequence of disk models which have the same Tcg{r/rg) 
distribution. Since 



T,4ff cx MMr~^RR{r/rg) oc M-^M{r/rg)-^RH{r/rg) , 



(4) 
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the same Teff distribution is obtained for models with fixed 
M /M"^. If disks radiated as blackbodies, all the spectra 
of the sequence would be exactly the same, only vertically 
shifted in the absolute value of the emergent flux. Indeed, 
the spectra are similar in the long-wavelength (optical and 
IR) portion of the spectrum, but they differ appreciably 
in the UV and EUV spectral region. In particular, the Ly- 
man edges of hydrogen and He II change their appearance 
significantly. Note that the non-LTE and LTE models for 
lower black hole masses in this figure are very nearly the 
same. This is expected, as the average density p oc M/M"^ 
times some function of r/vg for radiation pressure domi- 
nated annuli. Hence for fixed M/M'^, p scales as M~^, 
implying that the lower black hole mass models in figure 
15 have higher average densities and should therefore be 
closer to LTE. 



(M^.Mdot) pairs with the same cos i - 0.1 




frequency (Hz) 

Fig. 16. — Integrated spectral energy distribution uLu (in erg 
s~^) for the pairs of (Mg, M) which have the same L/i^Edd = 0.15. 
The numerical values of Mg and M expressed in Mq yr~^, are 
the same; the upper curve corresponds to the pair (32,32), and the 
lower one to (1/8, 1/8). Full lines: non-LTE models; dotted lines: 
LTE models. The models with the highest Mg (with the highest 
IR, optical, and UV flux) have the largest He II Lyman jump, and 
consequently the lowest flux for u > 2 X 10^®. 

Finally, in figure ^ we present a sequence of models 
with the same L/L^m- Since L/L^dd oc M/M, we chose 
a sequence where Mg and M (in Mq yr~^) have the same 
value; the corresponding L/Ledd = 0.15. Again, the shape 
of the spectrum is quite similar in the optical and IR re- 
gions, while there is a progressively larger portion of EUV 
radiation for lower black hole masses. The non-LTE effects 
are extreme for high-mass holes in the He II Lyman con- 
tinuum, for which the LTE models predict virtually zero 
flux. 

3.3.1. Effects of changing the viscosity parameter a. 

In figure ^ we compare the predicted spectra for disk 
models with Mg =1 for two values of the viscosity param- 
eter a. Although the spectra exhibit some differences, the 
effect of different a is very small in the optical and UV 
region; the only appreciable differences are found in the 
He II Lyman continuum. This result is very encouraging 
because it shows that the effect of the ad hoc viscosity 
parameter a is rather small, and therefore the predicted 



spectra are not influenced significantly by this inherent 
uncertainty. Similar conclusions were reached in Paper II, 
albeit for a few representative annuli. The sense of the 
net effect is that larger a leads to a lower density and 
thus larger departures from LTE, which cause a somewhat 
higher flux in the He II Lyman continuum. 
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Fig. 17. — A comparison of predicted spectral energy distribution 
for models computed with a = 0.01 (solid line), and with a = 0.1 
(dashed lines), for models with Mg = 1/2 and various values of M. 
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Fig. 18. — A comparison of predicted spectral energy distribu- 
tions for the sequence of models for a maximum rotating Kerr black 
hole (solid line), and for a Schwarzschild black hole (dashed lines), 
for models with Mg = 1 and various values of L/L^^^. The values 
of the mass accretion rates are 2, 1, 1/2, 1/4, etc., for the Kerr hole, 
and the corresponding values for the Schwarzschild hole are 5.613 
times larger. 

Finally, we present several representative spectra for the 
case of a Schwarzschild black hole. Again, the full set of 
spectra is available upon request. 

In figure ^ we compare the predicted spectra for disk 
models with Mg =1, for the maximum rotating Kerr and 
Schwarzschild black holes. The total luminosity of the 
corresponding pairs of models is equal; the mass accre- 
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tion rate is therefore a factor of 5.613 higher for the 
Schwarzschild case to adjust for the different efficiencies. 
The Kerr spectra tend to extend to higher frequencies for 
several reasons all arising from the fact that their disks 
extend in to smaller radii. As a result, the maximum ef- 
fective temperature found in the disk is higher, and the 
relativistic effects strengthening the high frequency spec- 
trum away from the axis are also greater. 

This completes our presentation of the overall integrated 
spectra of our grid of models. We now address some of 
their observational implications. 

3.4. Comparison with Other Models 

We first compare the spectrum of an LTE model with 
that computed by Sincell & Krolik (1998). Figure |l9| shows 
a comparison between one of their spectra and ours, com- 
puted for the same parameters. The agreement is satisfac- 
tory, as the spectra differ by at most 20% near the peak. 
The differences that do exist may be due to a number of 
factors, ranging from technical numerical contrasts in the 
methods to the detailed physical assumption made in both 
papers. 



that local blackbodies also produce a significantly larger 
flux in the UV and optical region. At IR wavelengths both 
models coincide because both use local blackbody flux for 
the cool annuli. 
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Fig. 19. — Comparison of an LTE disk model with Sincell & 
Krolik (1998). Solid line is our model; dotted line is theirs. Model 
parameters are Mg = 0.27, a/M = 0, i/iEdd = 0-3, Tout = lOOrg, 
o = 0.1, and the disk is viewed face-on. 

In figure ^ we show a second comparison, this time to 
a Laor & Netzer (1989; Laor 1990) model, kindly sup- 
plied by A. Laor. We choose a model with Mg = 1, 
M = 1 Mq yr-\ (i.e., L/L^^d = 0.15), a = 0.01, and two 
values of cosi: 0.8 and 0.2. We have integrated the local 
disk spectra out to a cutoff radius ront/fg — 217.8 to agree 
with Laor's value. The predicted spectra are generally 
similar, although there are several interesting differencies. 
Our models produce larger flux in the immediate vicinity 
of the He II Lyman edge (due to non-LTE effects leading 
to an emission edge for the hottest annuli), but a lower 
flux for the highest frequencies (likely because Laor & Net- 
zer take into account the effects of self-irradiation of the 
disk). Our models produce lower flux in the UV and opti- 
cal regions, which is a consequence of the different vertical 
structure and of non-LTE effects. Although the Laor & 
Netzer models do not simply assume local blackbody flux, 
figure is nevertheless quite indicative because it shows 



Fig. 20. — Comparison of a non-LTE disk model (solid curves) 
with a Laor & Netzer (1989; Laor 1990) model (dotted curves). 
Model parameters are Mg = 1, L/L^dd = 0.15, a/M = 0.998, 
rout = 2n.8rg, a = 0.01, and cosi = 0.8 and 0.2. 

3.5. Optical/ Ultraviolet Colors 

A common criticism of accretion disk models is that if 
they are to produce substantial ionizing photon flux, then 
they should have blue optical/ultraviolet colors. This is 
based on the long wavelength, low frequency behavior of 
a blackbody disk, which has F^, oc , with (3 — 1/3. We 
address the issue of ionizing photon flux in section 3.8 be- 
low, but here we wish to point out that our disks have 
quite red optical/ultraviolet colors. Figure |2^ shows the 
logarithmic spectral slope (3 as measured between 1450A 
and 5050A, where (3 is defined by the two corresponding 
frequencies by 

Fu2 V 



(5) 



Our disk models have colors near the median value (3 = 
—0.32 for bright quasars (Francis et al. 1991). Indeed, 
even disks with local blackbody emission have such red 
colors for these accretion rates and masses (Koratkar & 
Blaes 1999). This is because the temperatures are cool 
enough that the 1450-5050A spectra are not in fact in the 
long- wavelength limit. Note that our model spectra can be 
somewhat bluer than blackbody disks, but they are still 
sufficiently red to explain the colors observed in bright 
quasars. 

Figure |2^ shows some representative Kerr disk models 
viewed at i = 37° compared to the Francis et al. (1991) 
composite quasar spectrum. The shorter wavelength com- 
posite spectrum of Zheng et al. (1997), scaled to match 
the Francis et al. spectrum at 1285 A, is also shown. This 
latter composite is thought to be more trustworthy than 
the Francis et al. composite at wavelengths shorter than 
the Lya/NV line because of corrections for intervening ab- 
sorbers. The models were chosen to have the right color 
based on figure El], and then a least squares fit was done 
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to determine a single multiplicative normalization factor. 
The fit used supposedly line-free continuum windows of 
the Francis et al. (1991) composite spectrum, as defined 
by their figure 7: 1283-1289A, 1321-1329A, 1455-1475A, 
2196-2208A, 2236-2246A, 3024-3036A, 3928-3936A, 4035- 
4045A, 4150-4220A, and 5464-5476A. Note that the mod- 
els shown have M/M^ ranging from 1/16 to 1/2, with 
the lower mass models having the smaller values of this 
quantity. This is consistent with the behavior shown in 



figure 15: models with fixed M/M have slightly bluer 



optical/UV colors for decreasing black hole mass. 




Fig. 21. — Optical/ultraviolet spectral slope between 1450A and 
5050A for Kerr disk models with 37° viewing angle. Solid curves 
show our models with viscosity parameter a = 0.01, while the long 
dashed curves show local blackbody models. From top to bottom, 
the curves are for the nine increasing black hole masses of our grid: 
Mg = 1/8, 1/4, 1/2, 1,2,4,8, 16,32. 

These fits demonstrate that, while it is easy to recover 
the overall red color of the Francis et al. (1991) composite 
spectrum, explaining the shorter wavelength far ultravio- 
let emission seen in the Zheng et al. (1997) composite is 
more problematic. Two of the models shown in figure |2^ 
do in fact bracket the Zheng et al. composite, but they 
turn over at the shortest wavelengths shown in the figure. 
This might be a problem in view of the fact that an ex- 
trapolation of the Zheng et al. composite joins up with 
the ROSAT soft X-ray composite of Laor et al. (1997), 
implying that there is no cutoff. However, it is also con- 
ceivable that some other choice of model parameters (disk 
inclination, black hole spin) might work better. Notice 
also that Zheng et al. suggested that in order to explain 
the Lyman continuum flux one has to invoke the presence 
of a Comptonizing corona with temperature about 4 x 10® 
K with optical depth of the order of unity. In any case, 
composite spectra made from many sources may of course 
be unphysical, but these fits probably give some indica- 
tion of how our models will fare in explaining data from 



individual sources. 

It is also noteworthy that the low luminosity models 
shown in figure 22 have rather strong spectral features in 
the region of the hydrogen Lyman limit, which are gener- 
ally not observed in quasars. However, the models with 
the lowest luminosities are probably not representative of 
the quasars that make up the composite spectra. We now 
address this important issue of Lyman edges in AGN. 




Fig. 22. — Representative fits of some Kerr disk models with 
viscosity parameter a = 0.01 and cosi = 0.8 to the Francis et al. 
(1991) composite quasar spectrum (solid line). The shorter wave- 
length composite spectrum of Zheng et al. (1997) is also shown 
(solid line with more noise). From top to bottom, the dashed lines 
correspond to models with (Mq,M) = (8,16), (4,8), (2,1), (1,1/8), 
(0.5,1/64), (1/4,1/256), and (1/8,1/1024), respectively 

3.6. The Lyman Edge Region 

Quasars and active galactic nuclei are observed to have 
almost no intrinsic spectral features near the hydrogen Ly- 
man limit (e.g. Antonucci, Kinney, & Ford 1989; Koratkar, 
Kinney, & Bohlin 1992), and this has been a longstanding 
problem with accretion disk models (e.g. Krolik 1999a, 
Koratkar & Blaes 1999). To quantify the strength of Ly- 
man edge features in our model disk spectra, we calculate 
the relative change in flux at ±50A across the edge accord- 
ing to 



F 
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(6) 



Hence AF\/F\ is positive for an absorption edge and neg- 
ative for an emission edge. Figure |2^ shows the variation 
of this quantity with accretion luminosity and disk incli- 
nation angle for all our models around a Mg = 1 Kerr 
hole. Figure 2j shows the same thing but for a fixed in- 
clination angle of 37° and varying black hole masses. As 
was already apparent from the overall spectra shown in 
previous sections, substantial Lyman absorption edges are 
present in all our low luminosity disk models at modest, 
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near face-on viewing angles. However, high luminosity 
models have greatly reduced edges, particularly for the 
lower mass black holes. Because of their lower effective 
temperatures, higher mass black holes require higher Ed- 
dington ratios before the absorption edges are removed. 
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Fig. 23. — Measure of the hydrogen Lyman edge strength as 
a function of accretion luminosity, in units of the Eddington lumi- 
nosity, for a = 0.01 accretion disks around Kerr holes with mass 
Mg = 1. The curves are labeled with the value of cosi. 

It is important to emphasize that the absorption and 
emission edges can be smeared out by the varying Doppler 
shifts and gravitational redshifts in the accretion flow 
around the black hole when the viewing direction is at 
least somewhat off-axis. This is illustrated in figures 
and which show the actual spectral energy distribu- 
tions of some of our models in the Lyman limit region. 
The high luminosity models, which are probably most rel- 
evant for observed quasars, show very little in terms of 
sharp changes in flux, and our edge strength parameter 
shown in figures ^ and ^ really reflects an overall spec- 
tral slope, not an emission edge, in this wavelength region. 
These high luminosity models are still capable of explain- 
ing the observed red colors of quasars, as shown in figures 
2l| and |2l. 
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Fig. 24. — Same as figure but for fixed inclination angle of 
37° and varying Kerr black hole mass. 
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Fig. 25. — Spectral energy distributions in wavelength near the 
hydrogen Lyman limit for fixed Schwarzschild and Kerr black hole 
mass of Mg = 1 and viewing angle of 60° . From bottom to top, the 
curves represent accretion rates of , 2~*, . . . , 2", 2^ Mq yr~^, re- 
spectively, except in the Schwarzschild case where we have dropped 
the three lowest accretion rates. (The 2~^ Mq yr~^ Kerr curve is 
the same as theji^irve with the strongest Lyman absorption edge 
shown in figure p3.) Each curve has had its fiux multiplied by a 
different constantfactor in order to maximize clarity with all the 
curves on the same plot. All models shown here have a viscosity 
parameter a = 0.01. 

Reduction of the Lyman edge feature is caused both by 
relativity and by summing over emission and absorption 
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edges from the individual annuli at different radii. Only 
relativity (i.e. Doppler shifts) can smear out a flux discon- 
tinuity, however. We have tried integrating spectra with- 
out the relativistic transfer function for the M — 10^ Mq, 
M — 1 Mq yr^^, cosi = 0.8 case and found that without 
relativity, a substantial emission edge discontinuity exists 
in the integrated spectrum at 912 A. Relativity is not suf- 
ficient to smear out absorption edges in the low luminos- 
ity models, because the absorption edges are very strong in 
the spectra of all annuli in such models. Instead, the edges 
are simply blue shifted or red shifted away from 912 A. 



these models. 

In conclusion, the disk-integrated theoretical spectra for 
high Eddington ratio (i/LEdd) disks do not show signifi- 
cant features at the Lyman edge at 912 A, for both Kerr 
and Schwarzschild black holes. The only associated fea- 
ture is a change of the slope of the Lyman continuum, 
blue shifted by ~ 100 — 200 A, depending on the inclina- 
tion, mass, and to some extent on the black hole spin. It is 
likely that even this feature will be affected by additional 
physics, particularly metal line blanketing, which we will 
address in a future paper. 
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Fig. 26. — Same as figure ^ but for viewing angle of 37°. 

Sharp changes in slope are present in the low luminos- 
ity Kerr spectra near ~ 750A for cosi = 0.5, and in most 
of the spectra near ~ 840A for cosi = 0.8. In both cases 
these are at substantially shorter wavelengths than the Ly- 
man limit because of Doppler blue shifts, but they might 
nevertheless be observable in quasar spectra. From our 
models, we have calculated the maximum local change in 
logarithmic slope in the region 812-1012 A , and the results 
are illustrated in figure 27 for a viewing angle of 37°. We 
choose a sign convention such that a positive value of the 
slope change indicates a spectrum that becomes steeper 
as the frequency increases, i.e., positive slope change is as- 
sociated with absorption at the edge. All models display 
a similar behavior at this viewing angle, so we only dis- 
cuss the M — 10^ Mq case in detail (cf. figure ^ and 
the bold curve in figure At both low and high accre- 
tion rates, the maximum slope change occurs at ~ 840A 
for this black hole mass. As the accretion rate dimin- 
ishes, A((ilnFx/dln A) grows, refiecting an increasingly 
stronger smeared absorption edge. The reverse is true at 
high accretion rates. At intermediate accretion rates, the 
maximum slope discontinuity shifts to ~ 900A, where the 
spectra change from positive to negative slopes, reflecting 
a slight maximum in the flux around this wavelength for 



Fig. 27. — Maximum local change in logarithmic slope in 
812-1012 A region, as a function of accretion rate, for a = 
0.01 accretion disks viewed at 37° around Kerr holes of masses 
(2-3,2-2, ...,2*,2^)x 10^ Mq. The bold cuive is the M = 10^ Mq 
case illustrated in the bottom half of figure Kg. 

3.7. Polarization 

In addition to computing spectra, we have also calcu- 
lated the polarization in our complete grid of models, and 
this information is also available on request. Once again, 
no ad hoc assumptions are made here: the polarization is 
computed exactly in the radiative transfer calculation. In 
order to keep the parameter space as simple as possible, 
we do not include the effects of Faraday rotation by mag- 
netic fields in the photosphere of the accretion disk, which 
can be important in modifying the polarization signature 
(e.g. Agol, Blaes, & lonescu-Zanetti 1998). 

Figure shows the degree and position angle of the po- 
larization for various inclination angles for our a = 0.01, 
M = 1 M© yr~^ disk models around a M — 10^ Mq Kerr 
hole. These results are quite similar to those of Laor, Net- 
zer, & Piran (1990). In particular, the plane of polariza- 
tion is parallel to the disk plane at optical/UV frequencies, 
but rotates at higher frequencies due to general relativis- 
tic effects. Our predicted polarizations are higher than 
those of Laor et al. (1990), and our polarization generally 
dips redward and rises blueward of each continuum edge 
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(cf. especially the hydrogen Balmer edge and He II Lyman 
edge) . These differences are due largely to our more care- 
ful treatment of the vertical structure, the radiation field 
anisotropy, and the overall effects of absorption opacity. 

Polarization near the Lyman limit of hydrogen has pro- 
duced considerable recent interest due to the observation 
of steep rises in some quasars (Impey et al. 1995, Ko- 
ratkar et al. 1995, Koratkar et al. 1998). For illustration 
purposes, we show the degree of polarization at a view- 
ing angle of 60° predicted by our models for Mg = 1 black 
holes in figure |2^. This figure should be compared to figure 
25, which shows the corresponding spectra. Cooler disk 
models generally produce large polarizations, even larger 
than that for a pure electron scattering atmosphere (2.25% 
at this viewing angle, Chandrasekhar 1960). The reason 
for this is the enhanced anisotropy (limb darkening) of 
the radiation field due to vertical thermal source function 
gradients (Blaes & Agol 1996). However, steep rises in po- 
larization, and a steep rise in polarized flux as observed in 
PG 1630-^377 (Koratkar et al. 1995), are not produced by 
our models. This is partly due to the smearing and rota- 
tion of the plane of polarization by the relativistic transfer 
function (cf. Shields, Wobus, & Husfeld 1998). 




Fig. 28. — Polarization degree and position angle (6) for a = 0.01, 
M = 1 M0 yr-i disks around a M = 10^ Mq Kerr hole, for dif- 
ferent viewing angles. Solid curves correspond to our models, while 
dotted curves are those of Laor, Netzer, & Piran (1990). From top 
to bottom at high frequencies in the degree of polarization figure, 
the different curves correspond to cosi = 0.2, 0.5, and 0.8, respec- 
tively. This order is reversed (i.e. bottom to top) in the position 
angle figure. A position angle of 6 = 90° corresponds to the plane 
of polarization being parallel to the disk plane. 

The optical/UV polarization shown in figure Eq (in de- 
gree, position angle, and wavelength dependence) is gen- 
erally not observed in AGN and quasars (Berriman et al. 
1990, Antonucci et al. 1996). Once again, it is likely that 
the polarization of the radiation field will be affected by 
additional physics. In addition to Faraday rotation (which 



usually suppresses polarization), the additional absorption 
opacity from metal line blanketing in this region of the 
spectrum will probably reduce the polarization. Dust and 
electron scattering at larger distances from the continuum 
source can also modify the polarization signature (Kartje 
1995). 




Fig. 29. — Degree of polarization in the hydrogen Lyman limit 
region for disks around lO'' Mq Schwarzschild and Kerr black holes 
and a viewing angle of 60? . The models shown here are the same as 
those plotted in figure GSl with the bold curves in each plot repre- 
senting the highest and lowest accretion rates. (The lowest accretion 
rate gives the highest polarization around 800 A.) 

3.8. Ionizing Continua 

The radiation from accretion disks is often thought to 
supply most of the photoionizing continuum for the broad 
and narrow line regions of AGN. In view of this important 
application, we present the number of photons in the HI 
and Hell ionizing continua for disks with a range of ac- 
cretion rates and inclination angles around 10^ Mq Kerr 
holes in figures 30 and For comparison, we also show 
the ionizing continua for the corresponding models with 
local blackbody emission. For this particular black hole 
mass, our models generally predict somewhat fewer ioniz- 
ing photons in the hydrogen Lyman continuum than the 
corresponding blackbody disks. The exception is for near 
face-on disks at high luminosities. The reason is that at 
each annulus in the disk, our models generally have fewer 
low energy photons and more high energy photons com- 
pared to a blackbody at the same effective temperature 
(cf. figure 1^ and section 3.1 above). This can reduce the 
hydrogen Lyman continuum, while at the same time in- 
creasing the Hell Lyman continuum. Indeed, figure |3^ 
shows that, except at low luminosities, our disk models 
generally produce more Hell Lyman continuum photons 
than local blackbody disks. The reason for the dearth of 
photons at low luminosities compared to blackbody mod- 
els is the strong absorption edges present in these models. 
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0.4 0.5 0.6 
cos i 

Fig. 30. — Hydrogen Lyman continuum photon flux (times 4-Kd?, 
where d is the distance to the quasar) for a = 0.01 disks around a 
Mg = 1 Kerr black hole, as a function of viewing angle i. Solid and 
dashed curves correspond to our non-LTE models and local black- 
body disks, respectively. From top to bottom, the curves correspond 
to accretion rates of 2^, 2", . . . , 2~^, Mq yr"-*^, respectively. 
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Fig. 31. — Same as flgure ^ but for the Hell Lyman continuum. 
The lowest accretion rates are off the bottom of the figure. 

Note that the hydrogen Lyman continuum is Hmb dark- 
ened at high accretion rates, peaks at intermediate view- 
ing angles at moderate accretion rates, and is extremely 
limb brightened at low accretion rates (due to the high 
Doppler shifts overcoming the intrinsic absorption edges) . 
The Hell Lyman continuum is also limb brightened at all 
but the highest accretion rates, where it peaks at interme- 
diate viewing angles. These intrinsic anisotropics in the 



ionizing continuum produced by the disk may have impor- 
tant implications for photoionization models of the broad 
line region of AGN (cf. the much simpler anisotropy model 
of Netzer 1987) and for the statistics of AGN samples (e.g. 
Krolik & Voit 1998). 

4. CONCLUSIONS 

We have presented a grid of AGN disk models for a 
wide range of basic parameters, the black hole mass and 
the mass accretion rate, and for two values of the viscosity 
parameter a (0.01 and 0.1), and two values of the black 
hole angular momentum: maximum rotation Kerr black 
hole with a/AI — 0.998, and Schwarzschild black hole with 
a/M = 0. 

The basic aim of the present study was to construct a 
benchmark grid of models, based on simple, "classical" 
approximations, to which all our future, more elaborate 
models will be compared. The most important physical 
approximations that define the classical model are the fol- 
lowing: 

(i) The energy is generated by turbulent viscous dissi- 
pation, with the vertically-averaged viscosity described 
through the Shakura-Sunyaev parameter a. 

(ii) The vertical dependence of kinematic viscosity is de- 
scribed through a two-step power law in the column mass. 

(iii) Convection and conduction are neglected. 

(iv) No external irradiation or self-irradiation of the disk 
is considered. 

(v) Electron scattering is treated as coherent (Thom- 
son) scattering, i.e., the effects of Comptonization are ne- 
glected. 

(vi) Thermal opacity and emissivity of H and He only are 
included here, i.e., the effects of metals are neglected. 

(vii) Effects of line opacity are neglected. 

The underlying assumption, not listed here, is that the 
1-D approach is appropriate, i.e., that the disk may be 
described as a set of mutually non-interacting, concentric 
annuli. 

What is, however, treated exactly in this study, is the 
simultaneous solution of all the structural equations, with- 
out making any approximations concerning the behavior of 
the radiation intensity. Likewise, no a priori assumptions 
about atomic level populations (e.g. LTE) are made. The 
local electron temperature and density, mass density, and 
atomic level populations, are determined self-consistently 
with the radiation field. Once the local spectra of all an- 
nuli are computed, the spectrum of the full disk is found 
by integrating the emergent intensity over the disk surface 
using our relativistic transfer function code. 

We can summarize some of the overall spectral features 
of our benchmark grid as follows. Compared to multi- 
temperature blackbody accretion disk models, our spectra 
generally have lower fiuxes at low frequencies and higher 
fluxes at high frequencies. This difference is amplified fur- 
ther by relativistic effects, which are strongest for edge-on 
disks in Kerr spacetimes. Disks with different accretion 
rates around different mass black holes do not exhibit the 
same spectral energy distribution even if they have the 
same effective temperature distribution. Spectral slopes 
in the optical/UV region are significantly redder than the 
canonical F^, oc v^^^ low frequency accretion disk spec- 
trum. Non-LTE effects are important in all but the highest 
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density disks, enhancing the He II Lyman continuum and 
gencrahy reducing the strength of features near the HI Ly- 
man Hmit. HI Lyman edge discontinuities are only present 
in the cooler, low luminosity disk models. High Eddington 
ratio models exhibit no discontinuities, but they do show 
sharp changes in spectral slope, albeit at wavelengths sub- 
stantially blue shifted from 912 A. We stress again that 
by neglecting convection, the "cool" models {Tes < 9, 000 
K may be significantly altered; therefore the predicted 
optical and IR contimmm flux (and, in particular, the 
Balmer edge region) should be used with caution. Finally, 
our models show substantial wavelength dependent opti- 
cal/UV polarization which is parallel to the disk plane, 
a result which is likely to be modified by the effects of 
Faraday rotation and additional sources of opacity, both 
of which can suppress this polarization. 

In future papers of this series, we will systematically 
relax more and more of the classical assumptions listed 
above. Among these, relaxing assumptions (vi) are (vii) is 
straightforward, since our computer program TLUSDISK 
is fully capable of handling these situations. The only 
concern is that generating such models will require much 
more computer time than required for the present models. 
Also, one has to collect a large amount of atomic data, but 
we will profit enormously from already existing collections 
made for the purposes of modeling stellar atmospheres, or 
from data being included in current photoionization codes. 
Relaxing approximation (v) is less straightforward, but we 
have recently solved the problem and implemented Comp- 
tonization in TLUSDISK. Also, approximation (iv) may in 



principle be easily relaxed by adjusting the surface bound- 
ary conditions of each anmilus. 

However, relaxing assumptions (i) - (iii) is much more 
difficult. Here, we have to rely on detailed magnetohy- 
drodynamic simulations to guide us in how to describe 
convection, and how to choose the most appropriate pa- 
rameterization of viscosity. Nevertheless, even before such 
simulations are available, we can investigate phenomeno- 
logically the impact of a dissipation rate that varies with 
altitude within the disk. The existence of disk coronae 
suggest, for example, that the heating rate may increase 
with height, leading to a possible temperature inversion in 
the upper layers of disk atmospheres. 

Besides this purely theoretical motivation for improving 
disk models wc will also use the present grid of models to 
analyze a large volume of observed quasar spectra. Such 
a study will bring interesting results whether or not the 
models actually fit the data. If they do, this will be a 
strong argument in support of the accretion disk hypoth- 
esis, and of the adequacy of our theoretical description 
of accretion disks. If not, such a study will provide us 
with important clues as to which aspects of the theoret- 
ical description should be improved, and/or what other 
observational constraints will be needed in the future to 
settle these questions. 

We thank Ari Laor for providing us with his spectral 
models for comparison with those presented here. This 
work was supported in part by NASA grant NAGS- 7075. 



APPENDIX 

APPENDIX - DENSITY STRUCTURE OF THE DISK 

We consider here some details of the density structure in the case of a radiation pressure supported disk, without 
assuming that the gas pressure is totally negligible. 
We write the vertical hydrostatic equilibrium equation as (see Paper II for details of the formulation) 

^ + ^=-^^"' 

where Prad and Pgas are the radiation and gas pressure, respectively, p is the mass density, z is the vertical distance from 
the disk midplane, and Q = {GM/r^){C/B) [in the notation of Paper II; in the notation of Krolik (1999a), C/B = Rz{r)], 
where B, C and Rz{r) are the appropriate relativistic corrections. 

The radiation pressure gradient can be written (Paper II; Krolik 1999a) 

^ = -^F,.,, (2) 
dz c 

where Prad is the total (frequency-integrated) radiation flux, and xh is the flux-mean opacity. For most applications, the 
flux- mean opacity is well approximated by the Rosseland-mean opacity, kr. We consider here the case where the local 
opacity is fully dominated by electron scattering, in which case kr is constant and roughly equal to 0.34. (In fact, it is 
not exactly constant because it depends on the degree of ionization - see Paper II; however, we neglect this small effect 
here.) Further, we introduce (Krolik 1999a) 

= FU^) f{z) ^ asT^s /(^) , (3) 

where F°^^ is the total radiation flux at the surface, which is expressed through the effective temperature. 
We express the gas pressure through the sound speed, Cg, as Pgas = c^p. The sound speed is given by 

= T7 T, (4 

^ firriYi N -Tie 
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where k and mu are the Boltzmann constant and the mass of the hydrogen atom, respectively; fj, is the mean molecular 
weight (i.e., the mean mass of a heavy particle per hydrogen atom; in our case of a H-He atmosphere with a solar helium 
abundance /i = 1.4/1.1 = 1.27); N is the total particle number density, and ric is the electron density. The sound speed 
thus depends primarily on the temperature, and partly also on the degree of ionization [via the term N/(N — rie), which 
varies from 1 - for a neutral medium, to 2.1 - for a fully ionized solar-composition H-He plasma]. 
Substituting equations (||), (||), and (^) into (|^), we obtain 

1 del cl dp rj fi . 

where we introduce the radiation- pressure scale height, H^^ as 

which has the meaning of the disk height in the case of negligible gas pressure (see Paper II and Krolik 1999a). 

We shall consider two cases, (a) the gas pressure is completely negligible, and (b) the gas pressure is taken into account, 
although it is still smaller than the radiation pressure. 

In the former case, we take Pgas = 0, i.e. Cs = 0, and thus the l.h.s. of Eq. (|^) is zero. We are left with 

Hrf{z)-z = 0. (7) 

At the surface, /(zq) = 1 regardless of the dissipation law, and therefore zq = Hr, so the disk height is exactly equal to the 
radiation pressure scale height. However, the density cancels out exactly because both the gravity force and radiation force 
are linearly proportional to density (e.g., Krolik 1999a), so the density is undetermined by the hydrostatic equilibrium 
equation. In usual treatments of the radiation-pressure-supported disk, Eq. (|^) is assumed to hold everywhere in the 
disk, from which it follows that f{z) = z/Hr- In our approach, however, we assume that / is a known function of the 
column mass, m. In the case of constant kinematic viscosity (i.e., the dissipation rate proportional to density), we have 
(see Paper II), / = m/mQ, where mg is the column mass at the midplane, mg = S/2. This gives a linear relation between 
z and m, and since p — ~dm/dz, we obtain p{z) — const = po in the region of constant dissipation (99% of the total 
column mass of the disk in our models). In our numerical procedure, we in fact solve Eq. (|^) exactly, without assuming 
that the left-hand-side is negligible, and with / given as a known function of m (although not of z, because z is one of 
the state parameters to be solved as a function of m) . 

To be able to write down a simple analytic solution in the case of non-negligible gas pressure, we approximate the 
function f{z) as 

In other words, we assume that the radiation flux linearly increases with z until a certain height Hq, and then remains 
constant. Such a behavior is roughly observed in the numerical simulations. Since close to the midplane z/Hq = 1 — rn/mo, 
so that pq = —(dm/dz)z=o = 'm-o/Ho. Therefore, Hq = mo/po- Within the present analytical model, Hq, and thus po £^re 
quantities to be determined by the constraint of total column density (see below). We write 



where Hg has the meaning of the gas-pressure scale height corresponding to the midplane conditions (temperature and 
degree of ionization), and q{z) is a correction parameter that accounts for a dependence of Cs on z. The reason why Hg is 
called the gas-pressure scale height is that in the case of negligible radiation force (i.e. f{z) everywhere), the solution 
of Eq. (^ is given by 

p{z)-. pQexp[-{z/Hg)^], (10) 

where po = p{z = 0) is the density at the midplane. 

In the following, we neglect the dependence of sound speed on height, i.e. we assume q{z) = 1. One can derive 
appropriate analytical expressions even for a more realistic case at the expense of complicating the analytical formulas 
considerably, but the present approximation is satisfactory from the point of view of describing the basic physical picture. 
A similar analysis was already presented by Hubeny (1990). Equation (||) then reads, using Eq. (^, 

1 dp / Hr \ 2z 



p dz \ Hq J H^ 



1 — , for z < Ho (11) 



YTz=(^r-.)^ for.>Fo, (12) 
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which has the solution 



p{z) = po exp 



for z < H, 



, 



and 



The scale height is now determined 
after some algebra 











pocxp 




exp 


. Hg Hg 



for z > Hq . 

from the condition p{z)dz = mg. Substituting Eqs. dl^ ) and (|r 



where 



and 



(2/V7r) ho = y/ho/So erf (^\/ ho 6q^ + exp{~hr Sq) erfc((5o) , 



(13) 

(14) 

we obtain 
(15) 

(16) 
(17) 

where the error function erf(a;) is defined by erf(x) = {2/y/n) exp{—P) dt, and the complementary error function by 
erfc(a;) = (2/0F) exp{-t^) dt = l-erf(x). Equation (|l|) is a transcendental equation for Hq; however an approximate 
solution is found to be 

hoKihr + l/K. (18) 

If /if. 3> 1 (i.e., Hr ^ Hg), one obtains ho « hr, i.e., Hq « Hr, and thus p{z) « pQ for z < Hr', i.e., one recovers the 
solution for negligible gas pressure. If, however, Hg is not completely negligible with respect to Hr, we obtain Hq > Hr, 



ho — Hq/ Hg , hr — Hr / Hg 



Sn — hn — hr 



Hr/Ho 



1/2 



and density falls off exponentially with increasing height even in the midplane layer, with a scale height Hg/{1 
(see Eq. 

From the expression for the respective scale heights we see that Hr is weakly dependent on r {Hr oc D/C, so it depends 
on radial distance only through the relativistic corrections), while Hg oc (Tq/Q)^^^ where To is the temperature at the 



midplane. The latter scales as Tq cx mQ^^Tcff cx r and since Q oc r "^j we obtain finally H, 



3/8 



,,21/16 



i.e., it increases 



rapidly with radial distance. For the models displayed in figure I 
of complete dominance of radiation pressure), Hr — 2.63 x 10^ 
have Hq « Hr, which is indeed verified by the numerical model. 

to r/vg = 90, and we have Hr = 6.42 x 10^'^ cm, and Hg = 2.67 x 10^"^ cm, i.e., hr ~ 2.4. The correction 1/hr to ho is 
no longer negligible; we obtain /ip ~ 2.8, i.e., Hq « 7 x 10^^ cm (the exact value following from the numerical model is 
Ho — 7.7 X lO^'^ cm) and the density should show an exp[— (z/i?)^] decay with z with H ; 
roughly consistent with the numerical model. 



we have, e.g., at r/rg = 5 (which is well in the domain 
Hg = 5.52 X 10" cm, i.e., K « 48. We thus 
The last (coolest) model displayed there corresponds 
1013 



9.3 X 10^3 cm. This is indeed 
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